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Abstract 

We study the 3-state square-lattice Potts antiferromagnet at zero temper- 
ature by a Monte Carlo simulation using the Wang-Swendsen-Kotecky cluster 
algorithm, on lattices up to 1024 x 1024. We confirm the critical exponents 
predicted by Burton and Henley based on the height representation of this 
model. 
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1 Introduction 



Antiferromagnetic Potts models JT], |], H are much less well understood than their 
ferromagnetic counterparts. One reason for this is that the behavior depends strongly 
on the microscopic lattice structure, in contrast to the universality typically enjoyed 
by ferromagnets. As a result, many basic questions have to be investigated case- by- 
case: Is there a phase transition at finite temperature, and if so, of what order? What 
is the nature of the low-temperature phase? If there is a critical point, what are the 
critical exponents and the universality classes? Can these exponents be understood 
(for two-dimensional models) in terms of conformal field theory? 

One thing is known rigorously @, [§]: for q large enough (how large depends on the 
lattice in question), the antiferromagnetic g-state Potts model has a unique infinite- 
volume Gibbs measure and exponential decay of correlations at all temperatures, 
including zero temperature: the system is disordered as a result of the large ground- 
state entropyfj However, for smaller values of q, phase transitions can and do occur. 
One expects that for each lattice C there will be a value q c {C) such that 

(a) For q > q c (C) the model has exponential decay of correlations uniformly at all 
temperatures, including zero temperature. 

(b) For q = q c (C) the model has a critical point at zero temperature. 

(c) For q < q c {£) any behavior is possible. Often (though not always) the model 
has a phase transition at nonzero temperature, which may be of either first or 
second order. 

The problem, for each lattice, is to find q c (C) and to determine the precise behavior 
for each q < q c (C). 

For the common two-dimensional lattices, strong theoretical arguments^ — which, 
however, fall short of a rigorous proof — yield the following predictions for q c {C): 

{(3 + v / 5)/2 ps 2.618 ... for the hexagonal lattice 

3 for the square lattice /* ■ i \ 

3 for the Kagame lattice (L1) 

4 for the triangular lattice 

Monte Carlo simulations have confirmed numerically that the 3-state square-lattice 
model has a zero-temperature critical point ||, [7j , and that the 4-state square-lattice 
model || |7|] and the 3-state hexagonal-lattice model |8], |§ are non-critical uniformly 
down to zero temperature.^ 

1 This behavior has been proven for q > 4 on the hexagonal lattice, q > 6 on the Kagome lattice, 
q > 7 on the square lattice, and q > 11 on the triangular lattice ||. However, these bounds are 



presumably not sharp: see equation (1.1) below. 

2 Summarized in the introduction to fl. 

3 The Monte Carlo simulations of the 3-state hexagonal-lattice model reported in Ref. || give 
no evidence of any first-order phase transition as the temperature is varied from infinity to zero; 
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Two-dimensional models with zero-temperature critical points are of particular 
interest, as they can in most cases be mapped onto a "height" (or "interface" or "SOS- 



type" ) model |10|, O P], ff3, |T 



18, 19, 20, 21 



If this 



height model lies in its "rough" phase — a question that has to be investigated on 
a case-by-case basis — then its long-distance behavior is that of a massless Gaussian 
with some (a priori unknown) "stiffness" K > 0. The critical operators can then 
be identified via the height mapping, and the corresponding critical exponents can 
be predicted in terms of the single parameter K. In particular, if we know (by some 
other means) one of these exponents, then we can deduce the rest. 

Height representations thus give a means for recovering a sort of universality 
for some (but not all) antiferromagnetic Potts models and for understanding their 
critical behavior in terms of conformal field theory. All the nonuniversal details of 
the microscopic lattice structure are encoded in the height representation and in the 
stiffness parameter K. Given these, everything can be understood in terms of the 
universal behavior of massless Gaussian fields. 

The plan of this paper is as follows: In Section |2| we present briefly the general 
theory of height representations and then work out in detail the case of the 3-state 
square-lattice Potts antiferromagnet. Our presentation is based on the work of Henley 

supplemented by a few minor innovations of 



and collaborators |T9], ^ |2|, 0, 
our own. In the remainder of the paper, our goal is to test, by Monte Carlo simulation, 
the critical exponents predicted by Burton and Henley pB|| for the three relevant 
operators in the 3-state square-lattice Potts antiferromagnet at zero temperature. In 
Section [3] we describe our simulations, and in Section |] we analyze the data. 



2 Height Representations 



Many two-dimensional models with zero-temperature critical points can be mapped 
onto a "height" model: these include the triangular-lattice Ising antiferromagnet 



11, 12 



the triangular-lattice spin-S" Ising antiferromagnet | 2"3}| , the 3-state square- 

the 3-state Kagome-lattice Potts antiferro- 



lattice Potts antiferromagnet [10, 13, 26 



magnet ||T6| . F20| , the 4-state triangular-lattice Potts antiferromagnet ||27|| , the 4-state 
Potts antiferromagnet on the covering lattice of the square lattice [19, p0[ , a con- 
strained 4-state Potts antiferromagnet on the square lattice [f2B[, a special 6- vertex 
model [pOl, an d various dimer models 0, ffM [23, E5l and fully packed loop models 



[ |T7| , IS , |2"2"[ | . Here we shall explain briefly the basic principles underlying the con- 
struction of such mappings and their use to extract critical exponents. We shall then 
work out in detail the case of the 3-state square-lattice Potts antiferromagnet. 



this behavior is consistent with the theoretical prediction that the model has exponential decay 
of correlations uniformly down to zero temperature. However, these authors did not measure the 
correlation length or the staggered susceptibility, so no direct test of the non-criticality at zero 
temperature was made. Such a direct test is being made in Ref. 0. 
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2.1 General Theory 

The first step is to define a map assigning to each zero-temperature spin configura- 
tion {a(x)} a corresponding microscopic height configuration {h(x)}. This height rule 
is usually defined by local increments, i.e. one prescribes the change Ah = h(y) —h(x) 
in going from a site i to a neighboring site y in terms of the spin variables cr(x) and 
cr(y). For such a rule to be well-defined, one must verify that in all cases the net 
increment Ah around any closed loop is zero.[| The height variables h(x) lie in some 
discrete set TC C R D (for some suitable dimension D), which we call the height lattice. 

The next step is to identify the so-called ideal states: these are ground-state 
configurations (or families of configurations) of the original spin model whose corre- 
sponding height configurations are macroscopically "flat" (i.e. have zero net slope) 
and which maximize the entropy density (in the sense of maximizing the number of 
ground states that can be obtained from the ideal states by local modifications of the 
spins). We label each ideal state by its average height h G R D , and we define the 
ideal-state lattice X C R D to be the set of all average heights of ideal states. The 
equivalence lattice 

£ = {a ER D : a + I = 1} (2.1) 

is the subgroup of R D summarizing the underlying periodicity of X. 

We now guess that, in typical configurations of the spin model, the lattice is 
subdivided into reasonably large domains in which the spin configuration closely 
resembles one of the ideal states. It follows that typical configurations of the height 
model are given by domains in which the height h(x) exhibits small fluctuations 
around one of the values in the ideal-state lattice. We therefore expect that a suitably 
defined coarse-grained height variable h(x) will take values in or near the ideal-state 
lattice X, except at boundaries between domains. The long-wavelength behavior of 
the height model is thus postulated to be controlled by an effective coarse-grained 
Hamiltonian of the form 



H = \ dfx 



D 



vElv^l 2 + Mock(Mx)) 
1 i=i 



(2.2) 



where we have made explicit the components of the macroscopic height h = (hi, Ji2, ■ ■ ■ , 
The gradient term in ( |2.2|) takes into account the entropy of small fluctuations around 
the ideal states; the second term is the so-called locking potential, which favors the 
heights to take their values in X. We then expect that there exists some constant 
K r such that for K < K r (resp. K > K r ) the locking potential is irrelevant (resp. 
relevant) in the renormalizat ion-group sense. Thus, if K < K r our surface model is 
"rough" and its long-wavelength behavior can be described by a massless Gaussian 
model [28], with D components: 

(Mx) - h(y)] [h 3 (x) - h 3 (y)}) « l og \ x _ y \ (2.3) 



4 In fact, this property usually holds for free boundary conditions but not for periodic boundary 
conditions. We shall hereafter neglect this latter subtlety, by imagining that we are working always 
in infinite volume. 
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for \x — y\ ^> 1; in this case, the original zero-temperature spin system is critical. 
If, on the other hand, K > K r , then the surface model is in its "smooth" phase, 
exhibiting long-range order 

(h(x)) = h (2.4) 
and bounded fluctuations around this ordered state: 

[h(x) — h(y)}/ bounded as \x — y\ — > oo ; (2.5) 



correspondingly, the spin system is "locked" into small fluctuations around one of the 
ideal states. At K — K r the surface model undergoes a roughening transition.^] 

Let us note, finally, that a given ideal state can be represented by many different 
average heights h G X. More precisely, suppose that in some domain we have a 
particular ideal state X and that its average height is h G X. Now let us pass 
through various other domains of the lattice, coming back finally to a domain in 
which the ideal state is again X. We will find the average height in this latter domain 
to lie in the set h + TZ, where TZ is a particular subgroup of £ that we call the repeat 
lattice; we will also find, conversely, that whenever we enter a domain in which the 
average height lies in ho + TZ, that domain is in ideal state X. It follows that the ideal 
states are in one-to-one correspondence with the cosets X/TZ. 

The coarse-grained correlation functions of local operators in the spin language 
can be understood in terms of the correlation functions of local operators of the 
coarse-grained heights. The important point is that these latter operators should 
have the periodicity of the repeat lattice TZ. This means that the Fourier expansion 
of such an operator O, 

O(x) = £ G e iGrh{x) , (2.6) 

contains only wavevectors belonging to the reciprocal lattice of the repeat lattice, 

7Z° = {G G R D : G ■ a G 2ttZ for all a G 71} . (2.7) 

On the other hand, two wavevectors whose difference belongs to the reciprocal lattice 
of the equivalence lattice, 

£° = {G G R D : G ■ a E 2ttZ for all a G £} , (2.8) 

give rise to vertex operators exp[iG ■ h(x)} having identical long-distance behavior. 
The vertex operators of the height model are thus in one-to-one correspondence with 
the cosets 71° /£" ~ (£/7Z)°. 

Now, provided that the height model is in the rough phase (K < K r ), the corre- 
lation functions of the vertex operators exp[iG ■ h(x)] are given by 



(e^(°> e^*)) = 5 G , G ,exp -^([h(0) - h{x)f) 



6g,g< \x\- g2 '^ . 

(2.9) 



5 Note that the only alternatives for the spin model are criticality and long-range order. Thus, if 
there exists a height representation, the original spin model cannot be disordered at zero temperature. 
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It follows that the critical behavior of the operator O will be given by the most 
relevant vertex operator exp[iG ■ h(x)] appearing in its Fourier expansion with a 
nonzero coefficient: 

(O(O)O(ar)*) ~ (2.10) 

where 

G / v 

?7o = min . (2.11) 

(The scaling dimension is xo = rjo/2, and the operator is relevant in case the 
renormalization-group eigenvalue d — xq = 2 — xq is > 0.) This formula implies 
that we can write all the critical exponents in terms of a single parameter K. If one 
exponent is known, then all of them are. 

In particular, the locking potential V[ oc k has the periodicity of the ideal-state lattice 
X; its Fourier expansion ( |2.6|) has contributions only from wavevectors G belonging to 
8°. Let ago be the length of the smallest nonzero vector in 8°. Now, the roughening 
transition occurs exactly where the locking potential is marginal, i.e. where ??y lock = 4. 
It follows that 

* = (2 ' 12) 

If K < K r , the locking potential is irrelevant, with scaling dimension 

xy lock = 2K r /K > 2 . (2.13) 

It induces corrections to scaling ~ L 2 ~ xv ioek : where L is a suitable length scale. 

In addition to vertex operators exp[iG-h(x)], there is another type of local operator 
that makes sense in the massless Gaussian model: powers of gradients of h.^ In 
particular, the operator (Vh) 2n has scaling dimension 

x (V ft )2 n = 2n (2.14) 

and hence ?7(v/i) 2 « = 4n. It follows that all these operators are irrelevant, except the 
operator (Vh) 2 , which is marginal. Since these operators respect the lattice symme- 
tries, they can appear in the effective Hamiltonian and thereby induce corrections to 
scaling. The leading such operator is (Vh) 4 , with scaling dimension 2( V h) 4 — 4; it 
induces corrections ~ L 2 ~ X( - V ^ 4 = L~ 2 . 

Assuming that we have not overlooked any irrelevant operators that could ap- 
pear in the effective Hamiltonian, we conclude that the leading corrections to scaling 
behave as L _A , with 



(x Vlock - 2, x m , -2) = min^ - 2, 2) . (2.15) 



mm 



We remark, finally, the height representation can also be applied to these models 
at nonzero temperature. In that case one must consider also the fugacity of defects: 
that is, of places where the zero-temperature constraints are violated |2(], p7| . Very 
often the defect fugacity is a relevant operator. 

6 The height h itself is ill-defined as a field in dimension d < 2, due to infrared divergences. But 
gradients V/i, VV/j, . . . are well-defined. 
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2.2 Three-State Square-Lattice Potts Antiferromagnet 

The height representation of the 3-state square-lattice Potts antiferromagnet at 



zero temperature is very simple [10, 26]. Let the Potts spins a(x) take values in the 
set {0,1,2}. The microscopic height variables h(x) are then assigned as follows: At 
the origin we take h(0) = 0, 4, 2 (mod 6) according as <r(0) = 0, 1, 2; this ensures that 

h(0) = cr(0) (mod 3) (2.16a) 
h(0) = (mod 2) (2.16b) 

We then define the increment in height in going from a site x to a nearest neighbor 
V by 

h{x)-h{y) = a(x)-a(y) (mod 3) (2.17a) 
h(x) - h(y) = ±1 (2.17b) 

This is well-defined (in free boundary conditions) because the change Ah around any 
plaquette is zero{] It follows from (2.16) and (2.17) that 

h(x) = a(x) (mod 3) (2.18a) 
h(x) = x Y +x 2 (mod 2) (2.18b) 

for any site x = (x\, X2). In particular, the height h(x) is uniquely determined mod 6 
once we know the spin value cr(x) and the parity of x, and conversely. The height 
lattice TC is clearly equal to Z. 

There are six ideal states, given by 0/12 (spins on the even sublattice all equal to 0, 
spins on the odd sublattice chosen randomly between 1 and 2) and its permutations.^ 
In an ideal state, the height is constant on the ordered sublattice and fluctuates 
randomly ±1 around this level on the disordered sublattice. The average height of an 
ideal state is thus equal to its height on the ordered sublattice; it then follows from 
(2.18) that there is a one-to-one correspondence between ideal states and average 
heights mod 6 (see Figure |I])|] The ideal-state lattice X is thus also Z, as is the 
equivalence lattice £, while the repeat lattice 1Z is 6Z. The corresponding reciprocal 
lattices are £° = 2vrZ and TZ° = (tt/3)Z. 

There are three relevant operators (in the renormalizat ion-group sense) appearing 
in this model (see Table [I]) (T( 



M stagg (x) = (-ir +X2 a(x) (2.19) 

M u (x) = a(x) (2.20) 

Psta 99 (x) = \(-l) Xl+x > E (2^),^) - 1) (2.21) 

y nnn of x 



7 If four numbers ±1 add up to mod 3, they must necessarily be two +l's and two — l's, hence 
add up to 0. 

8 States like 0/1 are not ideal states because they do not have maximal entropy density. 

9 This fact was proven, in a different way, by Burton and Henley p6| , Appendix B.2]. 
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where we have represented the Potts spin at site x = (xi,X2) by a unit vector in the 
plane 

/ 27r 2tt \ 

a(x) = f cos — cr(x), sin — cr(x)J . (2.22) 

The first operator is the staggered magnetization; the staggering corresponds to a 
momentum k s4a9g = The second operator is the uniform magnetization. The 

third operator is a staggered sum over diagonal next-nearest-neighbor correlations 
(i.e. over y with \y — x\ = a/2); we call it the staggered polarization. In an ideal 
state, it takes the average value +1 (resp. —1) according as it is the even (resp. odd) 
sublattice that is ordered.0 

We can relate these observables directly to the microscopic height variables h(x) 
by exact identities. For the vertex operators with G = ±7r/3, ±27r/3, ±7r we have 

e ±iWM*) = M ^ agg {x)TiM^ agg {x) (2.23) 
e ±i(2n/3)h(x) = M^\x)±iM ( ?\x) (2.24) 

e ±inh(x) = ^_ 1 ^x 1 +« a (2.25) 



Here (|2.24 ) and (|2.25|) follow immediately from ( 2.18a|) and (2.18b), respectively, while 



( |2.23| ) follows by multiplying these two and taking the complex conjugate. Of course, 



the strictly local operator ( |2.25|) is trivial, but we can define a nontrivial almost-local 
operator with the same (G = 7r) long-distance behavior: 

e ±<W2)[fc(*)+k(y)] = (25^^) - 1) (2.26) 

for diagonal next-nearest-neighbor sites x, y^j It follows that G = tt corresponds to 
the staggered polarization. 

Remark. It is also of interest to define almost-local operators living on plaquettes. 
Let a; be a lattice site, and let = {(^i, ^2), (^1 + 1,^2), (a^i + 1,^2+1), (^1,^2 + 1)} 
be the plaquette whose lower- left corner is x. We then define the average height h(x) 
over that plaquette as 

h(x) = i £ %) . (2.27) 

yen(x) 

It is easy to see that h(x) takes values in Z U (Z + |): namely, it takes an integer 
(resp. half-integer) value if there are three (resp. two) distinct spin values a(y) on the 

10 Burton and Henley |2(| chose a slightly different definition of this operator: P stagg (x) = 
(l/4)(-l)^ £ S a{x)My) . 

y nnn of x 

11 Proof: Next-nearest-neighbor sites x,y always satisfy h(y) — h(x) = or ±2: the former case 
occurs when a{x) — cr(y), and the latter when a{x) ^ c(y). It follows that 



,±i(ir/2)[h(y)-h(x)) _ 2 r _ 1 

— ^ u a(x),a(y) L ■ 



Now multiply this by @.2E| ). 
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plaquette □(#). Indeed, the value of h(x) is uniquely determined mod 6 by the spin 
content of the plaquette: see Figure ||. Finally, for two adjacent plaquettes n(x) and 
□ (x'), we have 



Ah = h(x') - h(x) 



0,±|,±1 iih(x)eZ 
o,±| if^(x)eZ + | 



(2.28) 



The upshot of this construction is that, because h(x) takes half- integral as well as 
integral values, vertex operators exp[iGh(x)} and exp[iG'h(x)} are equivalent only if 
G = G' mod 4-7T, rather than mod 2tt as before; so we can define operators up to 
\G\ = 2n rather than only \G\ = 7r. (However, as will be seen below, all these "extra" 
operators are irrelevant.) We have 



±i(Tr/3)h(x) 



E M^ gg (y)TiM^ gg (y) 
yea(x) 



E Mlll gg {y)TiM2 gg {y) 

yea(x) 



E M£> gg (y)TiM£> gg (y) 

3 + (2^ — 3)5v( x ) >CT ( x ir}5 a (x') > a(a:''') 

(2.29) 



D ±i(2w/3)h(x) 



E MW{y)±iM®(y) 

y£0(x) 



E M i u\y)±iMZ , {y) 
yeu( x ) 



r(2), 



E M^{y)±iM^{y) 
yenjx) 

1 + $ar{x),ar{x")&a{x'),a{x'") 



;2.30) 



±iizh(x) 



1) 



X1+X2 



$<T(x),a(x") — 8ar(x'),a(x"') ± ^a{x),a{x")^cj(x l ),a(x"')A(a(x') 



a{x)) 
(2.31) 



^±i(4ir/3)h(x) 



e ±i(27r)h(x) e Ti(2n/3)h(x) 



E M«(y)T*M( 2 )(y) 
yen(x) 

1 — 35<j{x),(j(x'')fia(x>),a{x'") 



(2.32) 



^±i(5ir/3)h(x) 



e ±i(2n)h(x) e Ti(2n/3)h(x) 



E M^ fl (y)±iMj2k ff (y) 

3 — (2^ + 3)5 a (x),a(x")5a(x , )Mx'") 



(2.33) 



3 ±i(27r)h(a;) 



1 — ^cr(x),cr(x")Scr(x'),cr(x l ") 



(2.34) 



where we have labelled the sites around the plaquette D(x) as x,x',x",x"' in cyclic 
order, and in ( |2.31| ) we have used the shorthand A(n) = ±1 according as n = ±1 
mod 3. 



We can now read off the predictions for critical exponents. The staggered magneti- 
zation corresponds to G = tt/3 = a^o (this is the smallest nonzero vector in 7?°), hence 
VM.sta.gg = ir/(lSK). On the other hand, den Nijs et al. |10| and Park and Widom 
29| obtained the exact value T}M. stagg = 1/3 by means of a mapping to the 6- vertex 
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model. It follows that the height model corresponding to the 3-state square-lattice 
Potts antiferromagnet at zero temperature has stiffness K = tt/6. (In particular, we 
have K < K r = n/2, so the height model lies in its rough phase.) By the usual scaling 
law we obtain the susceptibility exponent (7 /u) stagg = 2 — r] Mstagg = 5/3. This value 
has been numerically verified by several authors |30], [], [JJ. 

The uniform magnetization corresponds to G = 2n/3 = 2an°- (In this model 
the ideal states have a nonzero net magnetization, which, however, is the same for 
A/BC and BC/A; the uniform magnetization is thus periodic on the ideal-state lattice 
with period 3.) It follows that r] Mu = 477 Mstoss = 4/3 and (7/1/) w = 2 - r/ Mll = 2/3 
|T0| . It is interesting that the uniform magnetization is predicted to have a divergent 
susceptibility in this antiferromagnetic model. We are not aware of any numerical 
test of this prediction in the literature. 

The staggered polarization corresponds to G = it = 3a^o. We have ^p sta99 = 
9VM atagg = 3 and hence {l/v)v stagg = 2- n Pstagg = -1 §. This means that the 
"susceptibility" for this operator does not diverge, but tends to a finite value with a 
power-law correction ~ L^ 1 (where L is the linear lattice size). This prediction has 
not, to our knowledge, been checked numerically in the literature. 

These are the only relevant vertex operators in the model. Indeed, a vertex 
operator exp[iG ■ h(x)] is relevant if and only if n = G 2 /(2ttK) < 4; or, writing 
G = na-jio with n integer, we need \n\ < V^nK/a-jio. The values K = ir/6 and 
a n° = tt/3 then imply that we must have \G\ < 2n/\/3, or \n\ < \/l2. 

The equivalence lattice has lattice spacing as = 1, so that the wavevector corre- 
sponding to the locking potential is G = 2n = 6a-^o and hence ^y lock = 36^M S i a9!? = 
12 > 4. So, V[ock is a (strongly) irrelevant operator. 

Remark. The foregoing predictions contain, at first glance, a serious paradox. 
The correlation functions of the microscopic staggered and uniform magnetizations, 

G stagg (x) = (M stagg (0)-M stagg (x)) (2.35) 

G u (x) = (M u (0) ■ M u (s)> (2.36) 

obviously satisfy 

Gst agg (x) = (-l) x ^G u (x) . (2.37) 

How, then, can G sta gg(x) decay at large \x\ like |x| -1 / 3 while G u (x) decays like |x| -4 / 3 ? 
The answer, presumably, is that the correlation functions contain both terms |10[ : 

G stagg (x) ~ Ixl" 1 / 3 + {-lf^\x\-^ + ... (2.38) 

G u (x) ~ (-l) x - 1+ * 2 |o:r 1/3 + |x|" 4/3 + ... (2.39) 

It is only when one passes to coarse-grained correlation functions, by smearing over 
several nearby lattice sites, that the oscillatory terms are replaced by much-more- 
rapidly decaying remnants, leaving 12 

G stagg (x) ~ \x\-^ + |xr 10 / 3 + ... (2.40) 

12 In order to get the maximum additional decay (namely, two powers of |ir|), it is necessary to 
smear over anmxn block with m and n both even. 
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G u (x) ~ |x|" 7/3 + |x|" 4/3 + ... (2.41) 

A similar cancellation of oscillatory terms occurs, of course, when one looks at the 
susceptibilities. 



3 Numerical Simulations 

In order to test all these predictions, we have carried out a Monte Carlo simulation 
of the 3-state square-lattice Potts antiferromagnet at zero temperature, on periodic 
L x L lattices with L ranging from 4 to 1024. We made our simulation using the 



Wang-Swendsen-Kotecky (WSK) cluster algorithm |3(J, [H]], which is ergodic at T = 
on any bipartite graph, and in particular on a periodic square lattice whenever the 
linear lattice size L is even p6|, 0.0 

For each lattice size, we made 10 6 measurements after discarding 10 5 iterations 
for equilibration. For L < 512 we performed a single long run starting from the 
ordered state 0/1. For L = 1024 we made two independent runs with different 
initial conditions, one starting in the ordered state 0/1 and the other starting in 
the ideal state 0/12 (each individual run was of total length 6 x 10 5 , with the first 
10 5 iterations discarded); there was no noticeable disagreement between the two sets 
of results. In units of the longest autocorrelation time Ti n t,p 2 t ( see below), our 
run length corresponds to ~ 1.3 x 10 5 Ti n t measurements, and our discard interval 
corresponds to ~ 1.3 x 10 4 r; n t iterations. This run length is sufficient to get a high- 
precision determination of both static and dynamic observables: we obtain errors of 
order < 0.2% for the static observables and of order < 2% for the dynamic ones.0 

Our program was written in FORTRAN. The runs for L < 512 were carried out 
on a Pentium 166 machine: each WSK iteration took approximately 5.7 L 2 /isec. 
The runs for L = 1024 were carried out on an IBM RS-6000/370 workstation, taking 
8.5 L 2 /isec per iteration. The total CPU time used in this project was approximately 
1 month on the former machine plus 4 months on the latter. 



13 By contrast, the WSK algorithm for q — 3 is known to be nonergodic on periodic 3m x 3n 
square lattices whenever m and n are relatively prime |32]| . Other cases are open questions. 

14 Our discard interval might seem to be much larger than necessary: 10 2 T; nt would usually be 
more than enough. However, there is always the danger that the longest autocorrelation time in 
the system is much larger than the longest autocorrelation time that one has measured, because 
one has failed to measure an observable having sufficiently strong overlap with the slowest mode. 
(Here is a minor example of this effect: the authors of Refs. ^j, [jj reported T- lnt < 5 because they 
failed to consider our slowest observable V s tagg, which has autocorrelation time T int ^ ss 8.) As 
an undoubtedly overly conservative precaution against the possible (but unlikely) existence of such 
a (vastly) slower mode, we decided to discard approximately 10% of the entire run. This discard 
amounts to reducing the accuracy on our final estimates by a mere 5%. 

Note also that while we have here performed our simulations only at zero temperature (/3 = oo), 
the authors of Refs. || Q employed a closely-spaced set of temperatures ranging from very high 
temperature ((3 = 2.0, £ « 5) to very low temperature ((3 = 6.0, £ « 20000) and found the 
autocorrelation times of M 2 stagg and the energy to be uniformly small. This constitutes further 
evidence against the existence of an undetected extremely slow mode. 
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The "zero-momentum" observables 

M s tagg = ^2 M ^agg(x) 

x 

M u = £M u (x) 

X 

'Pstagg = E Pstaqqjx) 

x 

all have mean zero. We have therefore measured their squares 

/ \ 2 



(3.1) 
(3.2) 
(3.3) 



stag 9 



Mt 



v 2 

stagg 



J2 M stagg(x) 



-E 

2V 



E M ^ 



x 



E 



x 

E ^<r x ,a 



2 



(3.4) 
(3.5) 



E P ' stagg{ x ) 



(xy) mm 



(3.6) 



as well as the "smallest-nonzero-momentum" observable associated to Nl stagg (x): 



stagg 



1 

2 



+ 



3 

- x 
2 2 



E(- 1 )" 1+:r2 e 27r " l/L ^ ! a 



+ 



E(-!r + * 2 e 27r " 2/ x,a 



(3.7) 



Here V = L 2 is the volume of the system, the sum ^ a is over the three possible values 
of the Potts spins, and the sum J2( X y) mm is over a U pairs of diagonal- next-nearest 
neighbors x, y (each pair taken only once). The staggered and uniform susceptibilities 
are given by 



Xstagg 



1 

V 



(M 2 sta J 



X U = y{Ml) 

and the "susceptibility" associated to the observable V sta gg is 

/(.Pstagg yV^stagg) ' 

Finally, the second-moment correlation length is defined by 

p [{.Xstagg / F stagg) ~ 1] ^ 

? ~ ' 2sin(vr/L) ' ~ ' 



(3.8) 
(3.9) 

(3.10) 
(3.11) 
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where ^ 

Fstagg = y stagg) ■ (3-12) 

The results of our simulations for the mean values of all these static observables are 
displayed in Table |2]. 

We have also measured the integrated autocorrelation time associated to each of 
the basic observables, using a self-consistent truncation window of width 6ri nt ]33 



Appendix C]. We find that the largest autocorrelation time (of the observables we 
measured) corresponds to Vg tagg , though all of them are roughly of the same order of 
magnitude (Table §). None of these autocorrelation times diverges as L grows; they 
tend to a constant. We have fitted the autocorrelation time for each observable to a 
constant (using methods to be described at the beginning of the next section). Our 
best fits are: 



lnt '' stagg 


= 7.552 ±0.052 




= 128, 


x 2 


= 5.43, 


3 DF) 


(3.13) 




= 4.921 ±0.027 




= 128, 


x 2 


= 1.87, 


3 DF) 


(3.14) 


T mt,M 2 t 

' stagg 


= 4.528 ±0.028 


(-^min 


= 256, 


x 2 


= 0.56, 


2 DF) 


(3.15) 


^"illt,^ stagg 


= 3.804 ±0.015 




= 32, 


x 2 


= 4.81, 


5 DF) 


(3.16) 



We conclude that the WSK algorithm for this model at T = has no critical slowing- 
down H [|J: r int < 8 uniformly in L. 



4 Data Analysis 

We perform all fits using the standard weighted least-squares method. As a pre- 
caution against corrections to scaling, we impose a lower cutoff L > L min on the data 
points admitted in the fit, and we study systematically the effects of varying L min on 
both the estimated parameters and the x 2 - I n general, our preferred fit corresponds 
to the smallest L min for which the goodness of fit is reasonable (e.g., the confidence 
level[^] is > 10-20%) and for which subsequent increases in L min do not cause the \ 2 
to drop vastly more than one unit per degree of freedom. 



4.1 Staggered Susceptibility 



The theoretically expected behavior of the staggered susceptibility at criticality 
(i.e., at zero temperature) is 



Xstagg 



(4.1) 



15 "Confidence level" is the probability that \ 2 would exceed the observed value, assuming that 
the underlying statistical model is correct. An unusually low confidence level (e.g., less than 5%) 
thus suggests that the underlying statistical model is incorrect — the most likely cause of which 
would be corrections to scaling. 
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with ( n f/v)stagg = 5/3; here A is a correction-to-scaling exponent and the dots indicate 
higher-order corrections to scaling. Based on the numerical results of Refs. [§, 0] , we 
do not expect large corrections to scaling on this observable. 

We tried first to extract the leading term in (|4.1|) by fitting our data to a simple 
power-law Ansatz Xstagg 

= ALW^^ss. This fit is reasonable already for L min = 32 
(X 2 = 4.34, 4 DF, level = 36%), but our preferred fit is L min = 128: 

i) 



= 1.66621 ± 0.00035 (4.2) 

V J stagg 

with x 2 = 1-31 (2 DF, confidence level = 52%). This result is only 1.5 standard 
deviations away from the expected value 5/3. 

We then considered the Ansatz Q4.1|), imposing the leading exponent (j/v) stagg — 
5/3 and trying various values for the first correction-to-scaling exponent A. We are 
able to find reasonably good fits already for L m i n = 4, provided we take A in the 
range 1.50 < A < 1.76. We therefore performed a three-parameter nonlinear weighted 
least-squares fit to simultaneously estimate A, B and A. Using L m i n = 4, we obtain 

A = 1.624 ±0.061 (4.3) 

with x 2 = 7.30 (6 DF, level = 29%). 

It is interesting to note that the exponent 5/3 is included in the interval (|4.3|) . 
If this is the true behavior, it means that the leading correction to pure power-law 
behavior in the staggered susceptibility is merely an additive constant: 

Xstagg = 0.87696(17)L 5 / 3 - 0.2820(30) (4.4) 

with x 2 = 7.78 (7 DF, level = 35%). Such a correction can be interpreted as a mere 
lattice artifact, not necessarily arising from any irrelevant operator of the continuum 
theory. 

4.2 Uniform Susceptibility 

The theoretically expected behavior for the uniform susceptibility is 

Xu = £ (7/iy) " [A + BL- A + ..] (4.5) 

with (7/1/) u = 2/3. The simple power-law Ansatz gives a decent fit only for L min = 
256, yielding 

fl\ = 0.6705 ±0.0022 (4.6) 

with x 2 — 0-35 (1 DF, level = 55%). This result is 1.75 standard deviations away 
from the theoretical prediction. 

The large deviations from pure power-law behavior for L < 256 can be explained 
as an effect of corrections to scaling. Indeed, if we consider the Ansatz ( |4.5|) with 
{l/v) u = 2/3 imposed and with just one correction-to-scaling term, we can obtain 
sensible fits even for L min = 4. But in this case, in contrast to the preceding one, 
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the range of acceptable A values is much narrower: 0.655 ^ A < 0.735. A three- 
parameter nonlinear weighted least-squares fit to A, B and A, with L min = 4, yields 

A = 0.695 ±0.013 (4.7) 

with x 2 — 1-96 (6 DF, level = 92%). In this case the value A = 2/3 is two standard 
deviations away from the above estimate, but the absolute discrepancy is small (less 
than 0.03) and can plausibly be explained as an effect of higher-order corrections to 
scaling. Indeed, the uniform susceptibility can be fitted well (with L min = 4) as a 
pure power law plus an additive constant: 

Xu = 0.54744(49)L 2/3 - 0.3486(16) . (4.8) 

with x 2 = 7.15 (7 DF, level = 41%). 

The results (|4.3D/(fOD and ([4.7p/(P~8"D, taken together, suggest that there are no 



irrelevant operators (having the symmetries of the Hamiltonian) with A < 5/3 and 
that the leading corrections to scaling in both Xstagg and Xu are lattice artifacts. This 
behavior is consistent with the prediction (|2.15|) that the leading irrelevant operator 
is (V/i) 4 , with A = 2. 



4.3 Staggered Polarization 

The finite-size-scaling behavior of Xp s ta gg is expected to be 

XP 3ta3g = XP sta J™) + BL~ A + ... (4.9) 

with A = 1. We tried first to ignore the correction-to-scaling term and fit the data 
to a constant. The fit is not very good: even for L min = 128 we have x 2 — 6.68 (3 
DF, level = 8%), with the estimate 

XP 3ta99 (oo) = 2.1736 ±0.0060, (4.10) 

and the confidence level gets slightly worse for L min = 256, 512. 

We next fit to (|4.9| ) with A = 1. For L min = 8 one already gets a fair (though not 
spectacular) fit: 

Xp sta Joo) = 2.1728 ±0.0054 (4.11) 

with x 2 — 9.63 (6 DF, level = 14%). However, the confidence level does not improve 
significantly for larger L min . 

We also tried fits to ( |4.9|) with various values A ^ 1. We were able to get reason- 
able fits for L min = 8, if we take 0.50 < A < 1.05. We then tried a three-parameter 
fit to obtain estimates for XPstaggi 00 )^ B and A. Our preferred fit corresponds to 

Xp sta Joo) = 2.160 ±0.011 (4.12a) 
A = 0.75 ±0.12 (4.12b) 

with x 2 = 6.39 (5 DF, level = 27%). The discrepancy between the above result 
and the predicted value A = 1 is only two standard deviations; it might be due to 
higher-order corrections. 
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4.4 Correlation Length 

Finally, we consider the scaling behavior of the second-moment correlation length, 
which is expected to be of the form 

f = U [x* + BL- A + . . .] (4.13) 

with p = 1. 

First, we tried to estimate the power p by a simple power-law fit. This gives a 
good result for L min = 128: 

p = 0.99875 ± 0.00069 (4.14) 

with x 2 — 0.48 (2 DF, level = 79%). This estimate is only 1.8 standard deviations 
away from the expected value p = 1, and the very small discrepancy (less than 0.0013) 
can be explained as an effect of corrections to scaling. 

If we look at Table |], we see that the ratio £/L increases from L = 4 to L = 8, 
decreases monotonically from L = 8 to L = 64, and then oscillates due to statistical 
noise for L > 64. Thus, if we want to study the L —>■ oo limit of this quantity 
without including correction-to-scaling terms, we expect to get a reasonable fit only 
for L min > 64. Indeed, if we fit our data to a constant x*, the first decent fit occurs 
for L min = 64, giving 

x* = 0.63546 ±0.00030 (4.15) 
with x 2 = 3.77 (4 DF, level = 44%). However, our preferred fit corresponds to 

-Lmin 512, 

x* = 0.63483 ± 0.00048 (4.16) 

with x 2 = 0.0051 (1 DF, level = 94%). 

On the other hand, if we want to study corrections to scaling, we must use at least 
some of the data with L < 64. The non-monotonic behavior for 4 < L < 64 indicates 
that, to obtain a reasonable fit over this whole interval, we would need at least two 
correction-to-scaling terms with amplitudes of opposite sign. An Ansatz with only 
one correction-to-scaling term could, at best, fit the data with L min > 8, and very 
likely not even that. 

Indeed, if we fit the data to the Ansatz ( J4.13Q with p — 1 and only one correction- 
to-scaling term ~ L~ A , we find that reasonably good fits are obtained for 0.25 < A < 
0.60 with L min = 8. (For L min = 4 we were unable to find any good fit, as expected.) 
We next tried a three-parameter fit to estimate x*, B and A. The first reasonably 
good fit corresponds again to L min = 8, and the estimates are 

x* = 0.63359 ±0.00132 (4.17a) 
A = 0.42 ±0.16 (4.17b) 

with x 2 = 9.35 (5 DF, level = 10%). However, a better fit is obtained with L min = 16, 
giving 

x* = 0.63479 ± 0.00066 (4.18a) 

A = 0.84 ±0.32 (4.18b) 
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with x 2 = 5.95 (4 DF, level = 20%). For L min > 32, we do not get any sensible result 
(A and B become very large, along with their error bars); this is due to the fact that 
most of these data correspond to the regime L > 64 where the corrections to scaling 
are submerged under the statistical noise. Let us remark that the value of x* given in 
( [4.18a| ) is only 1.6 standard deviations away from the one estimated by Ferreira and 
Sokal using extrapolation techniques at nonzero temperature: 

x* FS w 0.633888. (4.19) 

If we want to fit all the data (i.e. take L min = 4), we should introduce at least 
two correction-to-scaling terms. From the definition (|3 . 1 1|) , we expect two types of 
corrections to scaling for the correlation length: one of order L -5 / 3 coming from the 
numerator [cf. ( |4.4|) 1, and another of order L~ 2 coming from the subleading terms in 
the sine. Furthermore, we might also expect an effective constant-term "correction" 
of order L" 1 , analogously to what happened for the two susceptibilities. Thus, our 
next Ansatz would be 

j = x* + BL- 1 + CL- 5/3 . (4.20) 

If the coefficients B and C have different signs, the contribution of these two terms 
could be mimicked (in the range of monotonicity, L min > 8) by a single correction 
term with an exponent A e ff < 1, where A e ^j increases towards 1 as L min — > oo. 
Indeed, this scenario is in good agreement with our results (4.17b)/(4.18b). We 
therefore tried a three-parameter fit directly to the Ansatz ( [4.201) . Our preferred fit 
corresponds to L min = 4: 

x* = 0.63457 ±0.00033 (4.21a) 
B = 0.182 ±0.014 (4.21b) 
C = -0.521 ±0.035 (4.21c) 

with x 2 — 7.01 (6 DF, level = 32%). This result certainly does not prove that the 
Ansatz (|4.20|) is correct, since many other pairs of correction-to-scaling exponents 



could give an equally good fit; but it does display a satisfying agreement. 
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Operator 


G 


V 


'-f /u = 2 — i] 


Numerical Result 


M- Stag g 


±7T/3 


1/3 


5/3 


1.66621 ±0.00035 




±2tt/3 


4/3 


2/3 


0.6705 ±0.0022 


P stagg 


±7T 


3 


-1 


-0.75 ±0.12 



Table 1: Critical operators for the 3-state antiferromagnetic Potts model on the square 
lattice at zero temperature. The last column indicates the results from our Monte 
Carlo simulation (Section 



L 


Xstagg 


Xu 


%P-stagg 




4 


8.5576 ± 


0.0024 


1.0443 ±0.0010 


2.4486 ± 0.0076 


2.5127 ±0.0023 


8 


27.7671 ± 


0.0120 


1.8614 ±0.0028 


2.3888 ±0.0124 


5.1306 ±0.0055 


16 


88.8206 ± 


0.0441 


3.1537 ±0.0058 


2.2882 ±0.0128 


10.2497 ±0.0110 


32 


282.8085 ± 


0.1477 


5.2052 ±0.0105 


2.2464 ±0.0126 


20.4422 ±0.0219 


64 


897.3520 ± 


0.4776 


8.4362 ±0.0179 


2.2110 ±0.0124 


40.6662 ±0.0433 


128 


2851.2642 ± 


1.5484 


13.5899 ±0.0296 


2.1826 ±0.0120 


81.4388 ±0.0880 


256 


9056.7475 ± 


4.8888 


21.8141 ±0.0479 


2.1852 ±0.0122 


162.8235 ±0.1745 


512 


28731.6260 ± 15.5273 


34.7753 ± 0.0767 


2.1480 ±0.0116 


325.0478 ± 0.3449 


1024 


91167.3235 ±49.2042 


55.2604 ±0.1240 


2.1806 ±0.0120 


650.0264 ±0.6924 



Table 2: Mean values of the static observables for the 3-state square-lattice Potts 
antiferromagnet at zero temperature. 



20 



L 


T ^M 2 ataqq 


T int,Ml 


T "^ taqq 


int,^\ stagg 


4 


2.487 ±0.020 


1.525 ±0.010 


5.038 ±0.057 


3.330 ±0.030 


8 


4.073 ± 0.041 


3.240 ±0.029 


7.456 ±0.101 


4.032 ±0.041 


16 


4.417 ±0.046 


4.154 ±0.042 


7.845 ±0.109 


3.850 ±0.038 


32 


4.488 ± 0.047 


4.628 ±0.049 


7.765 ±0.107 


3.826 ±0.037 


64 


4.483 ± 0.047 


4.803 ±0.052 


7.863 ±0.110 


3.803 ±0.037 


128 


4.587 ±0.049 


4.903 ±0.054 


7.528 ±0.103 


3.855 ±0.038 


256 


4.558 ± 0.049 


4.908 ±0.054 


7.686 ±0.106 


3.816 ±0.037 


512 


4.515 ±0.048 


4.889 ±0.054 


7.377 ±0.100 


3.752 ±0.036 


1024 


4.512 ±0.048 


4.986 ±0.056 


7.638 ±0.105 


3.776 ± 0.037 



Table 3: Mean values of the dynamic observables for the 3-state square-lattice Potts 
antiferromagnet at zero temperature. 



L 


Z/L 


4 


0.62818 ±0.00057 


8 


0.64133 ±0.00069 


16 


0.64061 ±0.00069 


32 


0.63882 ±0.00068 


64 


0.63541 ±0.00068 


128 


0.63624 ±0.00069 


256 


0.63603 ±0.00068 


512 


0.63486 ±0.00067 


1024 


0.63479 ± 0.00068 



Table 4: Values of the ratio £/L for the 3-state square-lattice Potts antiferromagnet 
at zero temperature. 
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Figure 1: Ideal-state lattice X for the 3-state square-lattice Potts antiferro magnet. 
The symbols above the graph indicate the ideal states of the spin model: 0/12 means 
that the spins on the even sublattice are all equal to and that the spins on the odd 
sublattice are chosen randomly between the values 1 and 2. The numbers below the 
graph indicate the average height for the given ideal state; this height is determined 
modulo 6. 
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Figure 2: The average height h(x) on a plaquette is uniquely determined modulo 6 
by the spin content of that plaquette: 0/12 means, for example, that the two spins 
belonging to the even sublattice are both equal to 0, while the two spins belonging 
to the odd sublattice are 1 and 2. 
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